BCS - BEC crossover at T = 0: A Dynamical Mean Field Theory Approach 
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We study the T = crossover from the BCS superconductivity to Bose-Einstein condensation in 
the attractive Hubbard Model within dynamical mean field theory(DMFT) in order to examine the 
validity of Hartree-Fock-Bogoliubov (HFB) mean field theory, usually used to describe this crossover, 
and to explore physics beyond it. Quantum fluctuations are incorporated using iterated perturbation 
theory as the DMFT impurity solver. We find that these fluctuations lead to large quantitative 
effects in the intermediate coupling regime leading to a reduction of both the superconducting order 
parameter and the energy gap relative to the HFB results. A qualitative change is found in the 
single-electron spectral function, which now shows incoherent spectral weight for energies larger 
than three times the gap, in addition to the usual Bogoliubov quasiparticle peaks. 



1. INTRODUCTION 

The problem of the crossover from BCS superconduc- 
tivity to Bose-Einstein Condensation (BEC) of compos- 
ite bosons, where the superconducting coherence length 
(roughly the size of the fermion pair binding) is respec- 
tively much larger than or much smaller than the aver- 
age inter-fermion spacing, has been a problem of great 
interest from the very early stages of development of the 
theory of superconductivity. It was first addressed in 
the very early work of Eagles Q. In 1980 Leggett Q 
showed, using a variational approach, that at zero tem- 
perature the superconducting BCS ground state at weak 
coupling evolves smoothly into a Bose condensate state of 
tightly bound "molecules" at strong coupling. Nozieres 
and Schmitt-Rink Q extended the analysis to lattice 
models and to finite temperature and showed that the 
transition temperature T c between the normal and the 
superconducting state evolves continuously as a function 
of the magnitude of the attractive interaction between 
the fermions. The discovery of high T c superconductors, 
which are characterized by short coherence length com- 
parable to (but larger than) the inter-particle spacing, led 
to a resurgence of interest in the BCS-BEC crossover. A 
variety of interpolation schemes between weak and strong 
coupling developed using variational methods, functional 
integrals, and diagrammatic methods have been explored 
0,H,E1, and the existence of pseudo-gap anomalies in the 
normal state of a short coherence length superconductor 
has been established in two-dimensional systems 0, @ • 
Recently it has become possible to directly realize the 
BCS-BEC crossover in a dilute atomic gas of Fermions 
in a trap, by varying their two-body interaction (scatter- 
ing length) using a Feshbach resonance . 

In this paper we analyze the BCS-BEC crossover in 
the attractive Hubbard model using the dynamical mean 



field theory (DMFT) [T^, E3 approach. Our goal is to 
focus on the intermediate coupling regime U/t ~ 1 where 
one has no obvious small parameter. Since the DMFT 
becomes exact in the limit of infinite dimensions [ToL ll ll| , 
we are, in a sense, using the inverse coordination number 
of the lattice as the small parameter. The attractive 
Hubbard model has been studied recently using DMFT, 
but primarily in the normal phase |l2l| to analyze pair 
formation above T c and related phenomena. We focus 
here on the superconducting phase at zero temperature, 
in part because the DMFT method has been much less 
explored in broken symmetry phases. 

The remainder of the paper is organized as follows. In 
Section 2, we define the model and review the Hartree- 
Fock-Bogoliubov (HFB) mean field theory. In Section 
3, we briefly summarize the DMFT approach in the su- 
perconducting (SC) state and then describe the specific 
implementation of DMFT which we use, namely the iter- 
ated perturbation theory (IPT), in Section 4. In Section 
5 we present our results for the chemical potential, en- 
ergy gap, SC order parameter, density of states, spectral 
function, occupation probability and superfluid stiffness. 
We discuss how each of these evolves from the weak cou- 
pling BCS limit to the strong couping BEC limit, and 
to what extent the quantum fluctuations included in the 
DMFT implemented using IPT modify the results rela- 
tive to HFB mean field theory. 



2. MODEL 

We use the simplest lattice model which exhibits the 
BCS-BEC crossover, defined by the Hamiltonian 



H 



\u\J2 n *i n n 



(i) 



2 



The first term describes the kinetic energy of fermions 
with nearest neighbor hopping t, the on-site attractive 
interaction (— |Z7|) induces s-wave, singlet pairing and 
leads to a superconducting ground state for all b / 1, 
with the chemical potential /i determining the filling fac- 
tor n. (We will not study the system with n = 1 for 
which superconducting and charge-density wave orders 
coexist). 

The simplest mean- field description of this system uses 
the Hartree-Fock-Bogoliubov (HFB) theory leading to 
the following self-consistent equations for the "gap" A 
and /i at a temperature T = l/(fcs/3): 



and 
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We use standard notation where e^- is the band disper- 
sion for the fermions and £7. = \/A 2 + £{ w ith = 
e fc -/i- \U\n/2. As is well known BUS 13 at T = 
the solution of these equations leads to a BCS supercon- 
ductor in the weak coupling \U\/t <C 1 limit, to a BEC of 
hard core bosons in the opposite extreme \U\/t^> 1, and 
interpolates smoothly in between. However, the finite 
temperature solutions of these equations for \U\/t ^> 1 
yields a T c ~ |Z7|, not the BEC transition temperature 
scale expected to be of order t 2 /\U\, and the HFB ap- 
proach therefore does not constitute an interpolating ap- 
proximation at finite T Q. 

The DMFT is one of the simplest schemes that has the 
potential to overcome some of these limitations of simple 
HFB theory. As we discuss in the following sections, the 
lattice dependence of quantities that arise in the DMFT 
is not via the momentum k but only via the band disper- 
sion ek, and hence we can make the replacement 
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FIG. 1: The real part of the T-matrix T(w -> 0) for the two- 
body problem is plotted as a function of the attraction \U\/t 
between two fermions in an otherwise empty Bethe lattice 
of infinite connectivity. \U\/t = 2 is the threshold for the 
formation of bound state in vacuum. 



\U\/t where we expect the corrections to HFB at T = to 
be the most severe, and it will also emphasize the similar- 
ity between the continuum Fermi gases often studied the- 
oretically (and now experimentally) and the lattice model 
studied in this paper. The low energy scattering is de- 
scribed by the real part of the T-matrix ReT(w — > 0) for 
the two-body problem in vacuum, i.e., for two fermions 
in an otherwise empty lattice. This is the analog of the 
well known three-dimensional "scattering length" for the 
case of the Bethe lattice studied in this paper. As shown 
in Fig. ^ f° r \U\/t < 2 the attractive interaction is not 
sufficient to cause a two-body bound state in vacuum, 
and \U\/t = 2 is the threshold for bound state formation 
at which the scattering amplitude diverges [lij . We also 
note that at \U\/t = 2 the effective interaction diverges, 
i.e., one reaches the unitary limit, even though bare \U\ 
is in the intermediate coupling regime. We expect that 
the deviations from the HFB theory will be maximal in 
the vicinity of \U\/t — 2 where the system is efectively 
very strongly interacting. 



where p(e) is the (bare) band density of states(DOS). 
The implementation of the DMFT is often simplest on a 
Bethe lattice with a large coordination number z — ► oo, 
for which 
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if the bare hopping is normalized as 

t 



t 



(5) 
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We conclude this section by calculating the effective 
two-body interaction or the low energy scattering ampli- 
tude. This will give us a clear idea about the regime of 



3. DYNAMICAL MEAN FIELD THEORY 

To explore the intermediate coupling regime we use the 
dynamical mean field theory (DMFT) approach 
which reduces a lattice problem with many degrees of 
freedom to an effective single-site problem by "integrat- 
ing out" all the fermionic degrees of freedom except those 
at one site - the "impurity site" - and retaining the ef- 
fects of this only in the form of a self-consistently deter- 
mined bath with which the "impurity site" hybridizes. 
This retains non-trivial local quantum fluctuations miss- 
ing in conventional mean field theories and the descrip- 
tion can be shown to be exact in the limit of large di- 
mensionality. Since there are many excellent reviews of 
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DMFT we will only outline the elements of the technique 
in order to introduce our notation and to indicate the 
changes in the standard formalism necessitated by the 
presence of the superconducting long range order. 

To take the superconducting order (with singlet pair- 
ing) into account we use the Nambu formalism with the 
spinors = (cj^, c-fej.) and the matrix Green's function 



(5(fc,r) = -(T r *(fe,T)*t( fc) o)) 

G(k,r) F(k,r) 
Ft(fc,r) -G(-k,-7 



(7) 



where F(k,r) = — (T T Cfei-(r)c_fe;(0)) satisfies 
F(-k,-r) = F(k,r). We will denote all Nambu 
matrices by a 'hat' on top. In this formalism the 
interaction effects are described in terms of the self 
energy matrix 



S*(k, —iuj n 



S(k, iuj n ) 
-E*(k,iu n ) 



where uo n = (2n + l)ir/f3 are fermionic Matsubara fre- 
quencies. 

In the limit of infinite dimensions it can be shown that 
the self energy is purely local, i.e., is k independent (see 
ref. Hfl) . so that £ = T,(ioj n ) . Furthermore, the SC order 
parameter can be chosen to be real in a uniform system, 
which implies that S(iLu n ) = S* (—iuj n ) . Hence using the 
Dyson equation the full Green's function for the lattice 
can be written as 



G 1 (k,iw n ) 



ioj n + [i - e k 






lLU n — fl + Ck 



(9) 



iuj n + fj, - efe - E(icj„) -S(iu) n ) 

-S(iu n ) iu n - fjb + ek + T,(-iu n ) 



Thus, in DMFT the k-dependence in G~ 1 (k, iw„) enters 
only via the dispersion e k - 

The local self energy is itself obtained from an effective 
single site problem which can be regarded as arising from 
integrating out fermionic variables on all sites except one. 
The effective action for this single site problem within the 
DMFT approximation is given by 



S eff ~ 



dr / dT'^i^g- 1 ^ -t')#(t') 



-\U\ / drn i (r)n T (r). (10) 
Jo 

Here the host (Matrix) Green's function Q is not the non- 
interacting local Green's function, as it includes the ef- 
fects of the fermionic degrees at other sites which have 
been integrated out in the presence of interactions, i.e., it 



includes (local) self energy corrections at all these other 
sites, and needs to be determined by a triangle of self- 
consistency relations as described below. 

The first of these relations comes from the requirement 
that the "impurity" Green's function for the single site 
problem should be the same as the local Green's function 
of the lattice, so that 



G(iu! n ) = ^G(k,iuj n ). 



(11) 



This gives the diagonal and off-diagonal components of 
the impurity Green's function as 



G(iu> n ) -- 
and 

F(iuj n ) = S(iu> n ) 



de p(e) 



a 



(Ci -e){Q 2 -e) + S*{iuj n ) 
de p(e)- 



(12) 



1 



(Ci-e)(C2-e) + S 2 ^„) 
(13) 

Here (i = iuj n + p — T,(iuj n ) and (2 = iw„ £(— iw n )- 

For the case of the semi-circular DOS of eq. 10) the inte- 
grals can be evaluated in closed form as 



G(iuj v 



2Ci 



Xl + y/x\ - At 2 X 2 + \Jx\ - At 2 

(14) 



F(ico n ) = 



S(iuj n ) 
xi — x 2 



Xl + \Jx\- At 2 X 2 + yfx 



At 2 

(15) 

with 11,2 = C2 - Ci ± \/(C2 + Ci) 2 -45 2 (^„)/2 . 

The second relation comes from the Dyson equation 
connecting the full Green's function G at the impurity 
site, the host Green's function Q and the self-energy S, 
typically used in reverse, in the form 



Q (ioj n ) = G + Y,(iu n ). 



(16) 



The final relation comes from the solution for the self 
energy of the impurity problem defined by i|10|) . i.e., the 
determination of 



E(iw n ) = t,[g(iLu n )] 



(17) 



from a knowledge of the host Green's function. This is 
the task of the "impurity solver" , and is typically the 
hardest step in the triangle of self consistency. In this 
paper, we use iterated perturbation theory (suitably ex- 
tended to deal with the broken symmetry associated with 
the superconducting ground state as described in the fol- 
lowing section) as the impurity solver. 



4. ITERATED PERTURBATION THEORY 

We adapt the iterated perturbation theory (IPT), orig- 
inally developed for the paramagnetic phase of the re- 
pulsive Hubbard model |l0j,ll5|, to the SC phase of the 
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attractive Hubbard model. IPT is an approximate tech- 
nique which is much simpler than the more accurate 
but elaborate alternate methods such as quantum Monte 
Carlo [l6|. exact diagonalization [l7|. numerical renor- 
malization group |l8|. local moment approximation [l9j 
etc. IPT gives semi-analytical results which can be di- 
rectly and easily continued to the real frequency domain. 
It has been well studied in the context of the DMFT 
of the Mott transitions in the repulsive Hubbard model 
[lfl Hfil | where it gives results in complete qualitative 
agreements with the more accurate methods mentioned 
above, and only quantitative disagreement typically no 
more than 10-20 % in the transition temperatures and 
critical values of U/W; see, e.g., the comparison of re- 
sults obtained using different impurity solvers by Bulla 
et al. 18] . One can reasonably expect similar levels of 
qualitative and quantitative correctness in the present 
context in general. If qualitative changes are likely, this 
is commented on at appropriate places in the paper. 

In the form in which we use it here 20], IPT rests on 
the following ansatz for the self energy as a functional of 

G ■■ 



±IP T (U + ) = ±HFB + A£W(lU+). 



(18) 



Here Tjhfb is the Hartree-Fock-Bogoliubov (HFB) self 
energy as in eq. EH (see note [2l|), E^ 2 ) is the second 
order perturbation theory result (in powers of \U\) but 
calculated in terms of the Hartree- corrected host Green's 
function (see eq. I24|l , and A is to be determined as de- 
scribed below (see ea. 1281) . All the calculations we report 
and discuss in this paper are done at T = 0, and we work 
directly in real frequency oj + = lo + i0 + . 

The IPT ansatz is constructed so that [2(J it 

• reproduces the leading order terms for the self en- 
ergy in the weak coupling limit \U\/t << 1, 

• is exact in the atomic limit t/\U\ = 0, and 

• reproduces the leading order terms for the self en- 
ergy in the large u limit for all \U\/t, which ensures 
that some exact sum rules are satisfied. 

Thus IPT is expected to provide a reasonable interpolat- 
ing scheme between the weak and strong coupling limits. 
The HFB self energy is given by 



^hf = -\U\-t z -M x 



(19) 



Here t z and t x are Pauli matrices in Nambu space. The 
filling factor n = £ ct (4 c <t> and A = |C/|$ = |Z7|<CjC T ) 
with $ being the superconducting order parameter, are 
obtained from the full Green's function within IPT as 



A = -M / Im F{cj+)duj 

* J-oo 



(21) 



The diagonal and off-diagonal components of the sec- 
ond order self energy E^ 2 ) are given by 

E< 2 )(t) = -U 2 (&i(t)&2(-t)&aW - Fo\t)g 22 (-t)r (t)) 



(22) 



and 



S( 2 \t) = -U 2 (? {t)M-t)?o\t) - Gn(t)M-t)G22(t)) 

Here Gn, G22, 3q and Pq are components of the Hartree 
corrected Host Green's function matrix 

-1 , - , , _ , , v -1 



— E 



HFB 

(24) 

The subscript has been added in order to distinguish 
the components of the host green functions that arise 
in the specific context of the IPT approximation to the 
impurity problem. 

Each of the terms in Ij22(l and (|23|l is the product of 
three factors of the form 



H(t) = h 1 (t)h 2 (-t)h 3 (t) 



(25) 



where each hi is either £11, £22,^0 or ^o- Using the 
spectral representation for each Green's function we ob- 
tain for the Fourier transform, 



H(lu+) = - [ UdeiPiia) 

J — OO -_i 



N(ei,e 2 ,e 3 ) 



(26) 



where p l [e i 
mal factor 



- €1 + e 2 - e 3 
Im[hi(ef)]/ir and A^(ei, e 2 , 63) is a ther- 



N(e 1} e 2 , e 3 ) = f(ei)f(-e 2 )f(e 3 ) + f(- ei )f(e 2 )f(-e 3 ) 

(27) 

involving the Fermi function /(e) = 1/[1 + exp(/3e)] = 
l-/(-e). 

The matrix A in 118|) is fixed by demanding that 
E/p T (w + ) is "exact" in the large u> limit up to order 1/uj. 
As shown in Appendix A, we find A to be proportional 
to the identity matrix fo in Nambu space and given by 



.4 



U 2 n 



no 
2 



U 2 n 



TO 



(28) 

where fo is the identity matrix. Here no and Ao denote 
fictitious "filling factor" and "gap function" values eval- 
uated for the Hartree corrected Host Green's function, 
i.e., 



Im G{uj + )duj 



(20) 



riQ = —2/tt / hnQii(u) + )du) 



(29) 



5 



and 



A = -\U\/tt luLF a {uj)du;. 



(30) 



In the atomic limit, as discussed in Appendix B, we 
find that the second order self energy vanishes. Thus 
the IPT self energy for t/U = is simply the HFB self 
energy. We show in Appendix B that the HFB result is 
exact at zero temperature in the broken symmetry phase 
for t/U = 0, and thus our ansatz for the self energy is 
exact in the atomic limit. 
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5. RESULTS 

We have solved the DMFT equations within the IPT 
approximation as follows. For a given \U\/t and n, we 
start with a guess for the self energy and the 

chemical potential u . With this self energy as input, 
we compute the full local Green's function G(uj + ) us- 
ing analytically continued form of eqs. (|14|l and l|15l) at 
T = 0. Then we use the Dyson equation (|16|l to de- 
termine the host Green's function Q. Next we use the 
IPT ansatz I|18l24[l to determine the (new) self energy 
in terms of G, using new values of parameters A, Ao and 
no- Finally we obtain the new chemical potential by solv- 
ing the filling constraint equation (|2U[l using the Broyden 
method [22|. We then iterate the whole procedure until 
a self-consistent solution is reached, i.e., convergence in 
the self energy matrix is achieved. 

Within this self-consistency loop the evaluation of n 
and A using eq. (|2UI) and (|21[) (and similarly for no and 
Ao) involves integrals with singular integrands: the func- 
tions Im G(tt> + ) and Im F(w + ) have square root singu- 
larity at a gap edge uj = E g which is not a priori known. 
We fit these functions in a small neighborhood of the gap 
edge to the form K/ yJu> — E g , where the fits determine 
the gap in the spectrum E g . Then the singular part of 
integral near the gap edge is easily evaluated analytically 
and the part away from the singularity evaluated numer- 
ically using Gaussian quadrature. 

All of the results reported in this paper have been ob- 
tained at a fixed density of n = 0.5 (quarter filling) in or- 
der to avoid special features that arise at half-filling. (At 
half- filling, corresponding to n = 1, charge density wave 
order becomes degenerate with SC order and the Hamil- 
tonian has SU(2) symmetry, and is in fact isomorphic to 
the repulsive Hubbard model which has been well studied 
within DMFT.) We expect similar results for n ^ 1. 

Chemical potential: Fig. [21 shows the chemical po- 
tential u tuned to obtain n = 0.5 at T = 0. We see that 
it decreases monotonically as a function of \U\/t and the 
system becomes non-degenerate with increasing attrac- 
tion between the fermions. For \U\ > 3.5i the chemical 
potential goes below the bottom of the band. 



FIG. 2: The chemical potential n/t as a function of \U\/t for 
n = 0.5 and T = within IPT (filled circles) and HFB theory 
(full line). 
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FIG. 3: Single particle density of states (per unit energy per 
unit area) for n = 0.5 and T = within IPT for Bethe lattice 
of infinite connectivity. 

Density of states: The single particle density of 
states (DOS) N(w) 

N(u) = --Im G(w+) (31) 

is plotted in Fig. |3] for various values of \U\/t. We ob- 
serve a spectral gap (E g ) in the single particle DOS , 
which increases with U\/t as shown in Fig.0] where the 
gap as obtained within the HFB theory is also shown for 
comparison. 

For weak coupling the HFB spectral gap has the form 
t exp(— irt/2\U\), characteristic of BCS theory, while for 
large attraction it approaches to the binding energy of 
the composite bosons being proportional to \U\/2. The 
differences of the DMFT result for the energy gap from 
the simple HFB estimates will be discussed below. Near 
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FIG. 4: The spectral gap in the single particle density of 
states for n = 0.5 and T = as a function of \U\/t within 
IPT (filled circles) and HFB theory (full line). Note that the 
spectral gap within IPT is suppressed as compared to that 
obtained from HFB theory. 
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FIG. 5: The superconducting order parameter $ for n — 0.5 
and T = within IPT (filled circles) and HFB theory (full 
line). Note that the order parameter within IPT is suppressed 
as compared to that obtained from HFB theory. 

the gap edge the DOS has a square root singularity char- 
acteristic of a s-wave superconductor. But the DOS far 
from the gap edge does not simply look like the non- 
interacting semi-circular DOS of the Bethe lattice (as 
would be the case in weak coupling BCS theory). The 
structure at larger energy values comes from the ui de- 
pendence of the self energy, as we discuss below. 

Order parameter and energy gap: The supercon- 
ducting order-parameter is calculated using 



= -1/tt / ImF(w + )du 



and plotted in Fig. We see that, as expected, the 
quantum fluctuations included in DMFT suppress the 
order parameter below its HFB mean-field value. The 
effect of quantum fluctuations in the intermediate cou- 
pling regime can be seen more clearly in Fig. where we 
plot the fractional deviation of DMFT order parameter 



FIG. 6: 5$ = <&hfb — &ipt where <&hfb is the SC order 
parameter within HFB theory and $jpt is the same within 
IPT. SE g = (E 3 )hfb - {E g ) IPT where {E 9 ) H fb is the spec- 
tral gap within HFB theory and (E g )i P T is the same within 
IPT. 

and the energy gap from their corresponding HFB values. 
For small to intermediate values of the coupling U < t 
the quantitative differences are quite large with the HFB 
results being larger than the DMFT ones by more than 
100%. 

Spectral function: 

Another important quantity of interest is the one- 
particle spectral function 



A(e,u) = --ImG(e,u) + ) 

7T 



(32) 



where G is the "11" component of the Nambu Matrix 
Green's function for the lattice obtained by inverting eqn. 
|5J and is given by 



e, ijj = 



e — ji + £*(— up 



with 



(33) 



(34) 



[u + e - p + E*(-w+)] [w - e + n - E(w + )] - S 2 {u+) 

Since we are working within the DMFT framework, we 
have traded the k label for the energy label e. 

Quite generally, we expect that the spectral function 
will be of the form 

A(e, u) = Z+(e)S{w - E) + Z-{e)5(w + E) + A inc (e, u) 

(35) 

where Z± (e) are the coherent spectral weights in the Bo- 
goliubov quasiparticle/quasihole excitation poles at ener- 
gies ±E(e), and Ai nc is the incoherent part of the spectral 
function. We recall that in simple BCS-HFB mean field 
theory Z± = (1 ± £/E)/2 where £ = e - p - \U\n/2 and 
E = \J + A 2 and A,„ c = 0. In contrast, as shown in 
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FIG. 7: The spectral function A{e,uj) for n = 0.5, T = 
and \U\ = 2.5t as a function of uj for various values of £ = 
e — fi + E'(oj = E). Note that y4(e, uj) not only has coherent 
delta function peaks (which are shown by arrows with the 
corresponding weights) but also has broad incoherent parts 
which start appearing for uj > 3E g . For \U\ = 2.5t, E g — 
0.57i. 




FIG. 8: The imaginary part of self energy |E"(w)| for n — 0.5, 
T = and \U\/t = 2.5 as a function of w. Note that |£"(w)| 
becomes non-zero only for u> > 3E g 
0.57*. 



For |[/| = 2.5t, E g = 



Fig. [7| the DMFT result for A(e, uj) not only has sharp 
delta function peaks at ±E corresponding to the Bogoli- 
ubov excitations but also has a broad incoherent part. 
The weight in the coherent excitations Z + + Z_ < 1 and 
the deficit from unity is contained in A inc . All of this is 
a consequence of the frequency dependent self-energy as 
shown below. 

The imaginary part of the (diagonal) self energy is plot- 
ted in Fig.|Hl It vanishes at low energies (\ui\ < 3E g ) since 
there are no final states available for a scattering event. 
In this regime, from eqs. I|32|) and (|33|) it follows that 



A(e, w) = (u + e - (m + E(-w)) S (D(e, u)) . 



(36) 



with D(e, oj) = (uj + e — fi + £(— u>))(u) — e + (i — — 
S 2 (u>), since and S(u>) are now real. The quasipar- 



ticle excitation energies are the two symmetrically placed 
zeros at oj = ±2? of D which is an even function of uj. 
Thus the excitation energies are given by 



D{e,±E) = 0. 



(37) 



and the residues at the Bogoliubov quasiparticle poles 
are then reduced compared to their HFB values and are 
given by 



Z±{e) = [±E + e- (i + E(tE)] 



(38) 



The imaginary part of the (diagonal) self energy be- 
comes non-zero for uj > 3E g as shown in Fig. |SJ This 
can be seen to arise from the form of the second or- 
der self energy of eq. (J55Jl, because in a system with a 
gap, final states for scattering an injected particle off a 
particle-hole pair are available only if incident particle 
has uj > 3E g . We should note, however, that the 3E g 
value of the threshold is likely an artifact of DMFT/IPT 
which ignores collective excitations. It is well known 0] 
that in finite dimensions this model has a linearly dis- 
persing sound mode and scattering of a one-particle ex- 
citation off such a collective mode should lead to non-zero 
ImE(cij) above E g and not 3E g . 

In any case, within IPT, the structure of the self en- 
ergy leads to the incoherent spectral weight in A(e,uj) 
above three time the gap. The reduction of the coherent 
quasiparticle weight and the transfer of spectral weight 
to the incoherent part of the spetral function are thus 
features related to the uj dependent self energy and are 
missing in simple HFB mean field theory. 

Occupation probability: We next calculate the ana- 
log of the momentum distribution within the DMFT, 
namely the occupation probability n(e) of an energy level 
e given by 



n(e) = / A(e, uj)duj. 



(39) 



This is plotted in Fig. OH for various values of \U\/t. 

Within the HFB approximation, n(e) has the following 
simple form: 



n(e) = Z_(e) 



I - 



(40) 



It is easy to see that in the weak coupling limit n(e) 
looks like a slightly broadened Fermi function, dropping 
from 1 to zero over an energy scale of order A; its width 
hence increases monotonically with \U\/t. As one be- 
gins to form more and more tightly bound pairs, higher 
e states need to be involved in the pairing and eventually 
the system becomes non degenerate even at T = as 
already argued from the chemical potential. Note that 
n(e) within IPT is always less rounded than that within 
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FIG. 9: The occupation probability n(e) of an energy level e 
for n — 0.5, T — for various values of \U\/t obtained within 
IPT ( full curve) and HFB approximation (dashed curve). 
Top right panel shows n(e) for \U\ = 0.75t and top left panel 
shows n(e) for \U\ = l.Ot. Note e scale over which these 
are plotted are different. Bottom panel shows n(e) for \U\ = 
1.5i, 2.0t, 3.0i and i.Ot starting from bottom to top. 



HFB because quantum fluctuations reduce the gap in the 
single particle dos relative to HFB value. 

We find that the exact sum rule Jf° n(e)po(e)de = n/2 
is satisfied in our IPT calculation within the estimated 
numerical errors (0.4% — 2%). 

Superfluid stiffness: We can also estimate an upper 
bound on the superfluid stiffness D s which is the strength 
of the delta function in the real part of optical conduc- 
tivity : 



Reer(w) = D s 6(u>) + Kea reg (uj), 



(41) 



The Kubo formula for the superfluid stiffness 23] can be 
written as 



— = ~(JC X ) - ReA T (q x = 0, q y 

7T 



0,cj = 0) (42) 



where the kinetic energy — (tC x ) is the diamagnetic re- 
sponse of the system to the vector potential and the 
transverse current-current correlation function At is the 
paramagnetic response. It is easy to see that At > 0, so 
that D s < ir\(JC x )\. Thus the the kinetic energy gives an 
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FIG. 10: Upper bound on the superfluid stiffness D s for n = 
0.5, T — as a function of the coupling constant \U\/t within 
IPT (filled circles) and HFB theory (full line). 



upper bound to the superfluid stiffness, and in fact we 
may use it to provide a rough estimate of D s (although, 
as emphasized in ref . |24| , there is no reason to assume in 
general that for a lattice model D s is identical to 7r 1(^)1 
even though this equality holds within simple BCS-HFB 
theory.) In Fig. ^| we plot the superfluid stiffness D s as 
a function of the attractive interaction. We find that it 
is of order t in weak coupling, but decreases monotoni- 
cally with \U\/t reaching ~ t 2 /\U\ in the strong coupling- 
limit, which reflects the increasing effective mass of the 
hard core lattice bosons in the large \U\ limit. 



6. CONCLUSIONS 

In this paper we have studied the crossover from BCS 
superconductivity to BEC at T = in the attrac- 
tive Hubbard model using dynamical mean field the- 
ory (DMFT), implemented using the iterated perturba- 
tion theory (IPT) scheme. Our main goal was to ex- 
plore the DMFT approach in a broken symmetry state, 
which has received less attention than the paramag- 
netic phase, within a simple, easily implemented, semi- 
analytic scheme, and to see how the quantum fluctu- 
ations included in this framework modify the Hartree- 
Fock-Bogoliubov (HFB) mean field results. For the most 
part we found that HFB is qualitatively correct but over- 
estimates the SC order parameter and energy gap. In the 
intermediate coupling regime, the quantitative changes 
can be quite large. The frequency dependent self-energy 
of DMFT leads to the appearance of incoherent contribu- 
tions to the single particle spectral function at energies 
larger than three times the gap, and the consequent re- 
duction in the coherent spectral weight in the Bogoliubov 
quasiparticle/quasihole poles in the spectrum . 
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APPENDIX A : LARGE cu LIMIT 

In this Appendix we will first determine the large to 
limit of the self energy and then use it to fix the param- 
eter A in our IPT ansatz l)18[l. The uj — > oo limit can be 
obtained by a moment expansion of the Green's function 



G(uj + ) = 



1 



which follows from the spectral representation 

PG(t)de 



or 



(43) 



(44) 



and the definition of the n th moment of the density of 
states MC") = p G {u>)u) n dui. 

To evaluate the moments it is useful to go to a Hamilto- 
nian formulation for single site impurity problem, which 
requires introducing auxiliary degrees of freedom to de- 
scribe the "bath" . For the superconducting phase of the 
negative U Hubbard model, one possible choice for the 
impurity Hamiltonian is : 



H, 



imp 



■<£« 



ka 



ka 



v k (ctf ka + flc) +dJ2 (ftrf-ki + f-kifki 



(45) 

which describes the impurity c„ coupled to superconduct- 
ing bath of / fermions. Here Vk is the hybridization pa- 
rameter which allows fermions to hop between the bath 
and the impurity site and the D term represents s-wave 
pairing of the /'s. 

Using the spectral representation, the moments can be 
written in terms of commutators ([,]) and anticommuta- 
tors ({, }) involving the c's and the impurity Hamiltonian: 



(46) 



and 



< ^ = ({[c a ,H 1mp ],cl}) (47) 
= ([[[c a ,H imp ],H imp \c^) (48) 



where a, (3 —\, {. Explicit evaluation of these commuta- 
tors leads to the results 



M (1) = (-»-\U\n/2)T z +AT x 



M« = 



{2p\U\ + U 2 )n 



TO 



(49) 



(50) 



The Host Green's function Q is obtained from the im- 
purity Hamiltonian setting U — 0, and its large u> limit 
is 

g-'i^^w+ro + ^-Y.VkMrz (51) 



Using the Dyson equation, the large u> limit of self energy 
is then found to be 

«, ? U 2 n(l - n/2)/2 - A 2 , , 

E (uj+) = ^ HFB + f . 52 

U! + 

We must now find the large u> limit of IPT self energy, 
and compare it with the exact u> — > oo result l|52[l derived 
above. Begin with the diagonal component of E^ 2 ^(o;) 
given by 



£< 2 )(o 



= U" 



i=l 



9i (ei,e2,e3)iV(ei,e2,e3) 



iuf - ei + e 2 - e 3 



(53) 



where 



gi(ei,e 2 ,e 3 ) = Pii(ei)p22{e2)p22(e3)-pf(ei)p22(e 2 )pf(e 3 ) 

(54) 

Here pu(uj) = — 1/ir Im Qa{u) + ) with i = 1,2 and 
Pf{uj) = —1/tt Im Tq(lo + ). In the large uj limit it suf- 
fices to keep terms up to order 1/uj. We thus get 



E( 2 )(c+)^ 



^£ (1 _^)_ a 2 
2 v 2 ' 



(55) 



Next, consider the off-diagonal component of £^(a;) 



/oo «j 
f[de 



92(^1, e 2 ,e 3 )N(e 1 ,e2,e3) 



w+ — €1 + e 2 — £3 



(56) 



where 



02 (ei, e 2 , £3) = Pf(ei)p f (e 2 )pf(e 3 ) - pii{ex)pf(e 2 )p 22 {e z ) 

(57) 

It is easy to check that in the large ui limit vanishes 
up to order l/u>. 

Comparing the large w limits of the IPT ansatz 
Ejjfb + (uj) and the exact self energy Q52|). we find 



.4 



U 2 n 



' " 'l - H'l - A= 



TO 

(58) 



10 



APPENDIX B : ATOMIC LIMIT 



and 



In this Appendix first we first solve the attractive Hub- 
bard model exactly in the atomic limit t/U = and then 
show that our IPT ansatz for self energy is exact in this 
limit. In the atomic limit one can drop the hopping term 
in the Hubbard Hamiltonian (and also the hybridization 
term in the impurity Hamiltonian) so that 



H = -\U\n^ni - ^(n T + n x ) 



(59) 



The various sites decouple and so we have dropped the 
site label. The four states are |0), | f), I I) an d I Tl) with 
corresponding energies 0, — fi, —p and — 2/i — \U\. 

To study the broken symmetry phase we introduce a 
pairing field h, 



H = 



-\U\n ] n i - ^(n T + n{] - h{c^c\ + c x c T ) (60) 



and finally take the h — ► limit. The T = equations 
for the filling factor and order parameter are given by 



i + {h/xy 



and A 



-\U\h/\ 
1 + {h/Xf 



(61) 



where A = {-{2(1 + \U\) - ^J{2p + \U\) 2 + 4h 2 )/2 is the 
lowest eigenvalue of the Hamiltonian. We thus find that 
in the h — > limit we get the atomic limit solution: p = 
-\U\/2 for any filling n and A = \U\^n{2 - n)/2. 

Now consider the HFB equations in the atomic limit 
at T = 



and 



n = 2 



1 

W\ 



1 + 



p+\U\n/2 



yfc + WWW 2 + A 2 



2 v /( A1 +|[/|n/2) 2 + A 2 



(62) 



(63) 



The solution of these self consistent equations is p = 
-\U\/2 and A = | U\ yjn{2 - n)/2, which shows that HFB 
theory is exact in the atomic limit at zero temperature. 

Finally we will show that the second order self energy 
vanishes in the atomic limit. The Hartree corrected Host 
Green's function 



G (w+) = g- 1 {w+)-£ HFB 



reduces in the atomic limit to 
- -i 



g {w+) = w+T + {p + \U\n/2)r z + At x 



(64) 



(65) 



p f {w) = -A/\U\{6{w - \U\/2) 5{w + \U\/2)) (67) 

and substitute these in the expression for thcE( 2 V) and 
check that all the components of second order self energy 
matrix vanish in the atomic limit. 



It can be checked that the full Green's function in the 
atomic limit is identical to this, which means that & 2 ^ (w) 
vanishes. Alternatively one can calculate the density of 

states corresponding to Qo, which are given by 

pn{w) = {l-n/2)5{w-\U\/2)+{n/2)6{w+\U\/2) = p 22 {-uj 

(66) 
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